BranchedNcK <- function(PartTree,PartProbs,PartCosts,ParentIdx,ChildIdx,K) {
  
  # This function is designed to calculate the probability of N choose K 
  #   given a hierarchy within N.
  # We begin by identifying the second lowest level, then select those components
  #   from that second lowest level with children.
  # Now that each group of components at the lowest level is ID'd, collapse to the next level up.
  #   At each new level, ID the new components that are part of that group, 
  #   rinse and repeat and collapse to the next level up.
  # As the branches are collapsed, they will merge, eventually leading to the one
  #   true branch, so at each level we must recheck for branch commonality.
  #
  # All parts are considered MISSION CRITICAL - SM3 has no redundancy.
  # The probability is done based on replaceable/repairable parts. 
  # The rules are:
  #   1. The probability the component's Ghost Probability fails
  #      - The entire parent is replaced/repaired, sub-components are ignored.
  #      - This is simplest to begin with because it negates all other calculations
  #        down that branch.
  #   2. Combinatorics of Sub-Components (GP doesn't fail)
  #      - Count all the ways 1-K compoents can fail.
  #      - At the lowest level, this is simply N-choose-K
  #      - Any component that survives does not include it's children in N
  #   Note: The component probabilities themselves are only used to 
  #         compute 1-K probabilities of the parents. 
  #         This could one day be replaced by the derived probabilities.
  #
  # The relationship between component & sub-components is a little tricky, to say the least.
  # From one perspective the component probability includes all sub-component probabilities.
  #   However, given the MISSION CRITICAL requirement of every component, if any set of the 
  #   sub-components fail then the entire component fails.  
  # This begs the question: Should the sub-component probabilities add-up to the single overall 
  #   component probaility, or are they independent?  
  # If P(component) = f(P(sub-components)) = P(1-CG sub failures), 
  #   where, CG = total # sub-components 
  #     CG = C + G
  #     C = #subs worth RR (Repair vs. Replace)
  #     G = #subs not worth RR (i.e. Ghost Components)
  #     K <= C
  #   Then I can use these probabilities interchangeably. 
  #   If K = C :: P(component) = P(G)*P(1-C) :: i.e. everything must survive
  #   If K < C :: P(1-C) = P(1-K) * P((K+1)-C)
  #     Knowing P(1-K) is necessary to calculate the tipping point of 
  #       replacing the component-assembly vs. its sub-components
  #     Knowing P(G) is necessary to determine whether a component can be RR'd
  #       So when you use a component probability, you are naturally including the GP.
  #
  # Ghost Probabilities
  #   When you use P(1-K) you are naturally excluding GP to calculate the RR tipping point.
  #   This begs the question: Do I need to calculate the Ghost Probabilities?
  #     At each level, when I evaluate a sub-component failure, and especially the 
  #       first two levels, I am assuming I can RR that component.  
  #     This is consistent with the contractor's (Raytheon's) assessment that 
  #       the components listed in the RLA are worth RRing.
  #     GP and these hierarchical relationships are relevant when calculating the 
  #       probability that a component-assembly will fail.  In other words, rather than 
  #       using the component-assembly probability supplied in the RLA, which includes the GP,
  #       we use only the sub-component probabilities and calculate P(1-K). 
  #     Practically, when coding BranchedNcK.R, which evaluates two levels at a time,
  #       each component has a P(1-K) associated with it, which should always be used 
  #       in place of the RLA probability given because it does not include GP.
  #     Once I have P(1-K) I can estimate the GP, which must accompany each sub-component
  #       along with P(1-K).  
  #
  # Proability of Replacing a Component
  #   Bottom-Up
  #     Practically, P(GP) is treated as a single component, 
  #       which means RR requires P(1-k) * (1-P(GP)), 
  #       since when P(GP) we must Replace the Component
  #     P(ReplaceComponent) = P(RC) = P(GP) + sum(P(E>TP)), 
  #       where P(E>TP) = P(Events > $TippingPoint)
  #   Top-Down
  #     The GP for each component-assembly applies at that level, 
  #       but ALL it's parents' GP's must also survive for each sub-assembly failure combo.
  #       This means each event must be multiplied by the product of the parents' GPs.
  #     To do this, since each P(k) and the events within them are summed:
  #       P(RCb) = P(RC) *prod(RCp), where the sub b & p stand for Branch & Parents respectively.
  #   PREDICTION
  #     For the most indentured components, the requirement that the parents survive means 
  #       they will have lower probabilities.
  #     For the top components, the probability that a relativel small subset will cost more 
  #       than the entire component-assembly will also have lower and lower probabilities.
  #     There will be a sweet spot in the middle.
  #
  # Inputs
  #   PartTree,
  #   PartProbs,
  #   GhostProbs,
  #   ParentIdx,
  #   ChildIdx,
  #   NumBranchParts,
  #   K
  #
  # Output
  #   BNcK - list vector - Each BNcK[x] is a vector of probabilities with 
  #          a maximum of K and minimum of N sub-components for that part.
  #          BNcK is the top level of a two-tier list vector. The first tier contains each 
  #          level of the hierarchy, and the second tier contains the probility of 
  #          1-through-K failures in/down that branch.
  #
  # Written by: Sean M. Gonzalez, Willcor Inc., October 2012
  #___________________________________________________________________________________
  
  # Initialize Variables
  library(igraph)
#   source('~/Documents/Sean/Willcor/Aegis SM3/Sustainability/Risk Model/WillcorDSI-PPODS/ExtractKBranchCombos.R')
  pcird <- file("ProbCompIsReplaced.data", "a")
  
  maxLevel = max(PartTree) 
  N = length(PartTree)
  BP = list(maxLevel=maxLevel,N=N)
  
  if (maxLevel <2) { # Return that component probability
    warning("Number of levels is less than two (2)!!") 
    # Just do a traditional NchooseK and return it
    xx = logical(N)
    BNcK = sum(combn(N,K,function(kx) {
      xx[kx]=TRUE; 
      prod(PartProbs[xx]) *prod((1-PartProbs[!xx])) } ))
    return(BNcK)
  }
  
  # We begin at the bottom of the tree, where we know how to compute N-choose-K
  # from the component probabilities, and build our way up to the top.
  BNcK <- vector("list",N)
  Pg <- numeric(N)
  Ptp     <- Pg
  Prc_Bup <- Pg
  ConfPC  <- logical(N)
  # Populate the 'Tips' of the branches
#   branchTips = (PartTree == maxLevel)
#   tipIdx = which(branchTips)
  for (c in 1:N) { #tipIdx) {
    if ((length(ChildIdx[[c]]) ==1) && (ChildIdx[[c]] == 0)) {
      BNcK[[c]] = PartProbs[c]
      Pg[c]  = 0 #PartProbs[xx] # Probability Ghost - Pg
      Ptp[c] = PartProbs[c] # Probability Tipping Point - Ptp
      Prc_Bup[c] = Ptp[c]
      cat("Component #",c," - Probability =",BNcK[[c]],"\n",file=pcird)
    }
  }
  
  explored = 0
  minLevel = 1
  for (lvl in (maxLevel-1):minLevel) {
    cat("Evaluating Level - ",lvl,"\n\n",file=pcird)

    lvlLdx = (PartTree == lvl)
    lvlIdx = which(lvlLdx) # Index into PartTree, aka overall component index
    
    # Now I loop through this set of parents and calculate the probability of 
    #   K failures down that branch.
    # If there are no children for that component then the calculation is simple.
    for (p in 1:length(lvlIdx)) {
      xx = lvlIdx[p]
      children = ChildIdx[[xx]]
      # Ensure you are using the children ONLY at the next level down.
      clx = (PartTree[children] == (lvl+1))
      childLvlDown = children[clx]
      
      if (children==0) {
        numChildren = 0
        numChildBranches = 0
      }
      else { 
        numChildren = length(children) 
        numChildBranches = length(childLvlDown)
      }
            
      cat("\nEvaluating Component - ",xx," with children",children,"\n\n",file=pcird)
      
      # Calculate probabilities for 1-K
      # The Parent stands alone, but frankly shouldn't happen
      if (numChildren==0) {
#         stop("I thought I was precluding this...")
        BNcK[[xx]] = PartProbs[xx]
        cat("Component ",xx," BNcK =",PartProbs[xx],"\n",file=pcird)
      }
      # The sub-component stands alone.
      else if (numChildren==1) {
        BNcK[[xx]] = PartProbs[childLvlDown]
        cat("Component ",xx," BNcK =",BNcK[[xx]],"\n",file=pcird)
      } 
      # Each sub-component has depth
      else {
        # PartProbs changes depending on whether you're at the originating two levels, or 
        #   those components have their own sub-hierarchy to consider.
        # BNcK begins as empty, and because we are going up the hierarchy, when we call
        #   on the child probabilities, they are the latest ones to be populated.
        cat("\nFor component #",xx,"\n",file=pcird)
        cat("Calling ProbKSubCompFail - ",childLvlDown," - ",K,"\n",file=pcird)
        BNcK[[xx]] = ProbKSubCompFail(BNcK[childLvlDown],PartProbs[childLvlDown],K)
        cat("Component",xx," BNcK =",BNcK[[xx]],"\n",file=pcird)
      }
      
      # Can't calculate Ghost probabilities, child costs, etc. 
      # if there are no children
      if (numChildren==0 || all(children==0)) {
        Pg[xx]  = 0 #PartProbs[xx] # Probability Ghost - Pg
        Ptp[xx] = PartProbs[xx] # Probability Tipping Point - Ptp
        Prc_Bup[xx] = Ptp[xx]
      }
      # The sub-component stands alone.
      else if (numChildren==1) {
        # Calculate Ghost Probability using P(1-K)
        # P(component) = P(G)*P(1-C)  :: survival eq assumes P()'s are survival
        # P(component) = (1-P(G))*(1-P(1-C))
        #   P(G) = 1- [ P(component)/(1-P(1-C)) ] assumes P()'s are failures
        P1Ks = (1-BNcK[[xx]][min(numChildBranches,K)]) # Prob 1-K Survive -> P1Ks
        (1- ( (1-PartProbs[xx])/P1Ks) ) -> Pg[xx]
        # Here the probability of replacing the component is 
        # the Ghost Probability of the entire component because 
        # if the one sub-component fails, the only option is RR.
        Ptp[xx] = Pg[xx]
        Prc_Bup[xx] = Pg[xx]
        
#         cat("Single Child: Pg =",Pg[xx],"- xx =",xx,"\n")
        if (Pg[xx] <=0) { # Parent has greater reliability than the child!
          Pg[xx] = 0
          Ptp[xx] = min(PartProbs[xx:(xx+1)]) #I can go to xx+1 because there is only one child
          Prc_Bup[xx] = Ptp[xx]
          ConfPC[xx] = TRUE
#           tryCatch("Prc_Bup[xx] = Ptp[500000]",browser()) #Ptp[xx]
        }
      } 
      else { # Multiple children
        # Calculate Ghost Probability using P(1-K)
        # P(component) = P(G)*P(1-C)  :: survival eq assumes P()'s are survival
        # P(component) = (1-P(G))*(1-P(1-C))
        #   P(G) = 1- [ P(component)/(1-P(1-C)) ] assumes P()'s are failures
        P1Ks = (1-BNcK[[xx]][min(numChildBranches,K)])
        (1- ( (1-PartProbs[xx])/P1Ks) ) -> Pg[xx]
        
        # This is a check on the Ghost Prob because Raytheon did not design for this.
#         cat("\nMultiple Children: Pg =",Pg[xx],"- xx =",xx,"\n")
#         cat("Parent Prob =",PartProbs[xx],"- Branch Prob =",BNcK[[xx]],"\n")
#         cat("numChildBranches =",numChildBranches,"K =",K,"\n")
        if (Pg[xx] <=0) { 
          Pg[xx] = 0 
          ConfPC[xx] = TRUE
        }
        
        # Calculate where the RR Tipping Point is, then index it's probability (Ptp)
        # This could be done for each probability event, and that would be more precise
        SubCostAvg = median(PartCosts[children])
        SubCostStd = sd(    PartCosts[children])
        cat("Children =",children,": PartChildCosts =",PartCosts[children],": SubCostAvg =",SubCostAvg,"\n",file=pcird)
        # Compare cost of parent with sub-components
        Ktp = ceiling(PartCosts[xx]/SubCostAvg) # The K Tipping Point (Ktp)
        cat("Ktp =",Ktp,": PartCosts =",PartCosts[xx],"\n",file=pcird)
        
        thisDepth = length(BNcK[[xx]])
        if (Ktp > length(BNcK[[xx]])) { 
          Ktp = length(BNcK[[xx]]) 
          Ptp[xx] = 0
          blah = unique(diff(BNcK[[xx]]))
          blah = blah[blah>0]
          Prc_Bup[xx] = blah[length(blah)]
        } 
        else { 
          # Calculate Probability of Replacing Component
          # Assuming this component is the top of the hierarchy
          # P(ReplaceComponent) = P(RC) = P(GP) + sum(P(E>TP)), 
          #   where P(E>TP) = P(Events > $TippingPoint)
          blah = unique(diff(BNcK[[xx]]))
          blah = blah[blah>0]
          Ptp[xx] = blah[length(blah)]
#           Ptp[xx] = BNcK[[xx]][thisDepth] - BNcK[[xx]][max(Ktp-1,1)]
          Prc_Bup[xx] = Pg[xx] + Ptp[xx]
        }
      }
      explored = c(explored,xx)
      if (Prc_Bup[xx]<=0) {
        cat("Prc_Bup <=0!! Part #",xx,"with children:",children,"- Pg =",Pg[xx],"- Ptp =",Ptp[xx],"\n")
      }
      if (Pg[xx]<0) {
        cat("Pg <0!! Part #",xx,"with children:",children,"\n")
      }
      if (Pg[xx]>=1) {
        Pg[xx] = 0
        # This is the result of using average probabilities instead of iterating through all possibilities.
        # I would have to recode this so that I minimized variance by treating statistically
        # similar parts as one part.  How I would choose maximum variance IDK.
#         cat("Pg >=1!! Part #",xx,"with children:",children,"\n")
      }
      
      # LET'S MAKE A MOVIE!!
#       parentAssocGraph <- graph(ParentIdx$ParentAssocList, n=max(ParentIdx$ParentAssocList), directed=FALSE);
#       plot.igraph(parentAssocGraph,layout=layout.fruchterman.reingold,vertex.size=12) #,layout=layout.kamada.kawai)
      
    } # Parts within Levels
  } # Levels
  
  # Now that I have Probability of Replace from the Bottom Up (Prc_Bup), 
  #   let's add the top down.
  Prc = numeric(N)
  PKf = Prc
  cat("\nCalculating Total prob of comp replacement:","\n",file=pcird)
  for (c in 1:N) {
    px = ParentIdx[[c]]; px = px[px>2] # The use of 2 depends on variance of aggregated probabilities through the hierarchy
    PgParent = Prc_Bup[px]; PgParent = PgParent[PgParent>0]
    Prc[c] = Prc_Bup[c] *prod(1-PgParent)
    if (Prc[c]<=0) {
      cat("Prc =",Prc_Bup[c],"*",Pg[px],"\n")
    }
    
    lBNcK = length(BNcK[[c]])
    if (lBNcK ==0) {
      PKf[c] = 0
    } else { 
      PKf[c] = BNcK[[c]][lBNcK]
    }
    cat("\nComponent ",c,"Replacement Calcs:\n",file=pcird)
    cat("Prc = Prc_Bup *prod(1-Pg)\n",file=pcird)
    cat(Prc[c],"=",Prc_Bup[c],"*prod(1-",Pg[px],")","\n",file=pcird)
  }
  
  BP$BNcK    <- BNcK
  BP$PKf     <- PKf
  BP$Pg      <- Pg
  BP$Ptp     <- Ptp
  BP$Prc_Bup <- Prc_Bup
  BP$Prc     <- Prc
  BP$ConflictPC <- ConfPC
  BP$Explored <- explored
  
  return(BP)
}
#====================================================================================
ProbKSubCompFail <- function(BranchProbs,ChildProbs,K) {
  
  # The purpose of this function is to work at two levels only.  This function will
  # inevitably be called by other functions coordinating the branching.
  # For particularly large sets of children, I may employ some sort of average.
  #
  # Inputs
  #   BranchProbs - C x 1 asymmetrical list vector -
  #                 A list of the probability for each 1-through-K failures of branch 'c'
  #   ChildProbs - C x 1 vector - Probabilities that this component fails - from RLA list
  #
  # Outputs
  #   PartFailProb - K x 1 vector - Probabilities of failure for the parent component 
  #     for up to K potential failures.
  #
  # Written by: Sean M. Gonzalez, Willcor Inc., October 2012
  #__________________________________________________________________________________
  
  pcird <- file("ProbCompIsReplaced.data", "a")
  # Does not include the Ghost probability because this path requires that the 
  # Ghost failures do not happen.
  C = length(ChildProbs)
  cat("Number of Children for this component = ",C,"\n",file=pcird)
  
  # BranchDepths - C x 1 vector - Number of components that can fail for that branch
  BranchDepths = integer(C)
  cumDepth = 0
  for (c in 1:C) {
    BranchDepths[c] = length(BranchProbs[[c]])
    cumDepth = cumDepth + BranchDepths[c]
    cat("For c of C",c,C,": BranchProbs = ",BranchProbs[[c]],"\n",file=pcird)
  }
  
  # N-choose-K but with the hierarchical exclusions.
  # Whether the sub-component stands alone or has its own branch, ExtractKBranchCombos handles it.
  maxNumFailures = min(K,cumDepth)
  Prob1KFail = numeric(maxNumFailures)
  for (k in 1:max(maxNumFailures,1)) {
    xx = logical(C)
    
    # I NEED TO UPDATE THIS EXTRACTION FUNCTION TO GATHER THE PROBABILITIES AS
    # IT BUILDS TOWARDS K.  THE CURRENT METHOD REPEATS MANY OF THE STEPS AND IS
    # WOEFULLY INEFFICIENT!!
    # EACH K STANDS APART, AND SO WE WILL HAVE TO COMBINE THEM IN A LOOP THROUGH K
    # About the function itself, it is calculating the probability of 1-through-K failures
#     cat("Entering ExtractKBranchCombos.R \n")
    cat("\nEntering StarsAndBars.R - BranchDepths =",BranchDepths,": k =",k,": ChildProbs =",ChildProbs,"\n",file=pcird)
    cat("BranchProbs:",file=pcird)
    for (bp in 1:length(BranchProbs)) { 
      cat("\nBranchProbs",bp,"=",BranchProbs[[bp]],file=pcird)
    }
    cat("\n",file=pcird)
    
    if (k<=sum(BranchDepths)) { 
      PCNcK <- NchooseK_FlatBranchHybrid(BranchDepths,k,BranchProbs,ChildProbs)
  #     PCNcK <- StarsAndBars(BranchDepths,k,BranchProbs,ChildProbs)
  #     PCNcK <- ExtractKBranchCombos(BranchDepths,BranchProbs,ChildProbs,k)
    }
    else {
      cat("No more can be done with these PC relationships.\n",file=pcird)
      break
    } # No more to be done.
    cat("\nDone with branch perm K=",k,"\n",file=pcird)
    Prob1KFail[k] = PCNcK$Prob
#     PartFailProb[k] = sum(combn(C,k,function(kx) {xx[kx]=TRUE; prod(ChildProbs[xx,k])*prod((1-ChildProbs[!xx,k]))}))
    
    cat("For k =",k," of K total",maxNumFailures,"::\n",file=pcird)
    cat("PCNcK$NumEvents =",PCNcK$NumEvents,"\n",file=pcird)
    cat("PCNcK$Prob =",PCNcK$Prob,"\n",file=pcird)
    cat("Prob1KFail =",Prob1KFail[k],"\n",file=pcird)
    if (k != 1) {
      Prob1KFail[k] = Prob1KFail[k] + Prob1KFail[k-1]*(1-Prob1KFail[k])}
    cat("Prob1KFail = ",Prob1KFail[k],"\n",file=pcird)
    if (is.na(Prob1KFail[k]) || Prob1KFail[k] ==0) { stop() }
  }
  return(Prob1KFail)
}
